rm(list=ls())
setwd("C:/Users/bschneer/Documents/GOV1/GOV 2001/Replication/CEM")
library(cem)
library(mvtnorm)
library(modeest)
library(foreign)

files<-c("sexdisc2","sexhar2")
n<-0
#Determine Circuits to run models for
circuits<-c("1","2","3","4","5","6","7","8","9","10","11","DC")
circ_combn<-c()
for (i in 1:12) circ_combn[[i]]<-combn(circuits,i)

holder<-as.data.frame(matrix(0,nrow=4095,ncol=9))
names(holder)<-c("N Circuits","Circuits","ATT", "SE (ATT)","Beta","SE","L1 (Pre-CEM)","L1 (Post-CEM)","T in Sample After Match")

load("genderjudging.RData")

match<-get(files[1])[,c("judge.vote","treatment","circuit","dec.year","jcs","year.birth","minority.judge","lower.dir","jud.experience","abs.jcs")]
match["circuit"]<-as.factor(match[,"circuit"])

dec.year_cut<-c(1995, 1996.17, 1997.33, 1998.5, 1999.67, 2000.83, 2002)
jcs_cut<-c(-0.79, -0.33, 0.12, 0.58)
year.birth_cut<-c(1905, 1930, 1955)
minority.judge_cut<-c(0,1)
lower.dir_cut<-c(0,.5)
jud.experience_cut<-c(0,1)
abs.jcs_cut<-c(0, 0.53, 1.07)
circuit_gps<-list("1","2","3","4","5","6","7","8","9","10","11","DC")

cp<-list(dec.year=dec.year_cut,jcs=jcs_cut,minority.judge=minority.judge_cut,lower.dir=lower.dir_cut,jud.experience=jud.experience_cut,abs.jcs=abs.jcs_cut,year.birth=year.birth_cut)
gps<-list(circuit=circuit_gps)

for (i in 1:12){
for (j in 1:ncol(circ_combn[[i]])){
n<-n+1
new_circuits<-(circ_combn[[i]])[,j]
match1<-match[match$circuit %in% new_circuits,]

if (nrow(match1)>40){
imb1<-imbalance(match1$treatment,match1,drop=c("judge.vote","treatment"))
cem.match1<-cem(treatment="treatment",data=match1,cutpoints=cp,grouping=gps,drop=c("judge.vote"))
imb2<-imbalance(match1$treatment,match1,drop=c("judge.vote","treatment"),weights=cem.match1$w)

model<-glm(judge.vote ~ treatment,binomial(logit), data=match1,weights=cem.match1$w)

scale<-mean(fitted(model))*mean(1-fitted(model))

att<-scale*coefficients(model)[2]
att.se<-sqrt((scale^2)*summary(model)$cov.unscaled[2,2])
beta<-model$coefficients[2]
se<-sqrt(summary(model)$cov.unscaled[2,2])

holder[n,]<-cbind(i,paste((circ_combn[[i]])[,j],collapse=","),att,att.se,beta,se,imb1$L1$L1,imb2$L1$L1,cem.match1$tab[2,2]/cem.match1$tab[1,2])
}
#else break
}
}

output<-holder[order(holder$ATT),]

write.csv(output,file="sexdisc_combn.csv")